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Abstract 

A formulation of finite-rate ablation surface boundary conditions, including 
oxidation, nitridation, and sublimation of carbonaceous material with pyrolysis gas 
injection, has been developed based on surface species mass conservation. These surface 
boundary conditions are discretized and integrated with a Navier-Stokes solver. This 
numerical procedure can predict aerothermal heating, chemical species concentration, 
and carbonaceous material ablation rate over the heatshield surface of re-entry space 
vehicles. In this study, the gas-gas and gas-surface interactions are established for air 
flow over a carbon-phenolic heatshield. Two finite-rate gas-surface interaction models 
are considered in the present study. The first model is based on the work of Park, and the 
second model includes the kinetics suggested by Zhluktov and Abe. Nineteen gas phase 
chemical reactions and four gas-surface interactions are considered in the present model. 
There is a total of fourteen gas phase chemical species, including five species for air and 
nine species for ablation products. Three test cases are studied in this paper. The first case 
is a graphite test model in the arc-jet stream; the second is a light weight Phenolic 
Impregnated Carbon Ablator at the Stardust re-entry peak heating conditions, and the 
third is a fully dense carbon-phenolic heatshield at the peak heating point of a proposed 
Mars Sample Return Earth Entry Vehicle. Predictions based on both finite-rate gas- 
surface interaction models are compared with those obtained using B' tables, which were 
created based on the chemical equilibrium assumption. Stagnation point convective heat 
fluxes predicted using Park’s finite-rate model are far below those obtained from 
chemical equilibrium B’ tables and Zhluktov’s model. Recession predictions from 
Zhluktov’s model are generally lower than those obtained from Park’s model and 
chemical equilibrium B' tables. The effect of species mass diffusion on predicted 
ablation rate is also examined. 

Nomenclature 

B' = ml p e u e C m , dimensionless mass blowing rate 

C; = mass fraction for species i 

(C-i) = adatom i, O or N 

C m = Stanton number for mass transfer 

D = bifurcation diffusion coefficient, m 2 /s 

D, = diffusion coefficient for species i, m 2 /s 

E = total energy per unit volume, J/m 3 

F = nonlinear equation defined in Eq (24), or Pq / ^2mn~kT 
fi = diffusion factor of species i 
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h = Plank’s constant, J-s, or enthalpy, J/kg 
J = mass diffusion flux, kg/m 2 -s 
Kj = equilibrium constant 

K t = thermal conductivity of translation temperature, W/m-K 
K v = thermal conductivity of vibration temperature, W/m-K 

k = Boltzmann constant, J/K 

kf = forward reaction rate, defined in Eq (1 8) 

k r = backward reaction rate, defined in Eq ( 1 8) 

mi = mass of species i, kg 

m = mass flux, kg/m 2 -s 

M = molecular weight, kg/mole 

Nj = defined in Eq (12) 

p = pressure, N/m 2 

Pe = saturated vapor pressure, N/m 2 

q C onv = convective heat flux, W/m 2 

q v = heat flux due to species diffusion, W/m 2 

R = universal gas constant, J/kmol-K 

R b = base radius, m 

Rc = comer radius, m 

R n - nose radius, m 

r, = reaction rate, defined in Eq (14) 

S - recession rate, m/s 

T = temperature, K 

t = time, s 

u = fluid velocity, m/s 

v w = mass injection velocity, m/s 

x = Cartesian coordinate system, m 

Zj = bifurcation diffusion quantity of species i, defined in Eq. (4) 

a - surface absorptance 

P = efficiency of gas-surface interaction 

£ = surface emissivity 

£; = factor in i th heterogeneous reaction 

©° = free surface concentration 

©i = surface coverage concentration of species i 

X = blowing reduction parameter in Eq.(12) 

p = density, kg/m 3 

a = Stefan-Boltzmann constant, W/m 2 -K 4 
x = shear stress N/m 
/}, - defined in Eq. (4) 

q = general body fitted coordinate system normal to surface, m 
v = s pecies diffu sion flux, m/s 
^ = ^kT w !2mrii , m/s 
w = species source term in Eq.(8), kg/m 3 -s 
x = mole fraction 
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V = gradient, m' 1 
subscripts 
c = char 

E = chemical equilibrium 
e = boundary-layer edge 
g = pyrolysis gas 
i = species or direction component 
j = surface species, or direction component 
s = gas species or stagnation point 
v = virgin 
w = wall 


Navier-Stokes Solver 

The Navier-Stokes solver, GIANTS, is used to estimate the hypersonic 
aerothermal heating distribution over a blunt body. The governing equations used in the 
code, as developed by Lee, 1 may be characterized as representing a fiowfield in thermo- 
chemical nonequilibrium. The GIANTS code solves the time-dependent conservation 
equations of mass, momentum, and energy for the chemical and thermal non-equilibrium 
fiowfield. The species mass conservation equation is given by 2 


d Ps , d 
dt dxj 


( Ps u j)= ~ 



(Ps v sj)+ W s • 


( 1 ) 


the total momentum conservation is written as: 


d , d , , dry 

T-(P«i)+T— (pU|Wy)=-- — 
dt ox ,■ J ox 


and the total energy conservation is written as: 


( 2 ) 


8E _0_ 
dt dxj 


({E + p)u j ) = --^-(q‘j + q v j ) - ~ ( Ui Tij ) - Y.~ v ^ h 


dxj 


dxj 


s=\ dx j 


SJ'^S ' 


( 3 ) 


These equations are discretized using the implicit flux-split finite-volume method. The 
discretized equations may be solved by a block-tridiagonal matrix inversion using Gauss- 
Seidel line relaxation with alternating sweeps in the backward and forward directions. 2 " 4 
This technique has been shown to yield steady-state results efficiently. There is a total of 
fourteen gas phase species: C0 2 , CO, N 2 , 0 2 , NO, C 2 , C 3 , CN, H 2 , HCN, C, N, O and H. 
Five of the species are for air, and the rest are for ablation products. The gas phase 
chemical reactions implemented in the GIANTS code are from the work of Park. 5 The 
total number of reactions is nineteen, including nine dissociation reactions and ten 
exchange reactions, as follows: 

Dissociation reactions 
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(1) C0 2 +M ++ CO+O+M 

(2) CO + M <-> C + O + M 

(3) N 2 + M ++ N + N + M 

(4) 0 2 + M +-» O + O + M 

(5) NO + M <-+ N + O + M 

(6) C 2 + M ++ C + C + M 

(7) C 3 + M ++ C 2 +C + M 

(8) CN + M <-> C + N + M 

(9) H 2 + M ++ H + H + M 

Exchange reactions 

(10) NO + O <-> N + 0 2 

(11) N 2 +0 ++ NO + N 

(12) CO + Oe0 2 + C 

(13) C0 2 +0 +» 0 2 +CO 

(14) CO + C <-» C 2 +0 

(15) CO + N ++ CN + O 

(16) N 2 + C ++ CN + N 

(17) CN + O ++ NO + C 

(18) CN + C ++ C 2 +N 

(19) HCN + H ++ CN + H 2 

The mass diffusion models implemented in the GIANTS code include the constant Lewis 
number (Le), the constant Schmidt number (Sc), and the Bifurcation approximation. For 
the bifurcation model, the species mass diffusion rate is defined as 6 : 

J t =-pD^VZ i (4) 

Ml 

where 

7 Q ~ vQ - „ ,y Cjf 

%i ~ > Ml -lu , Hi =MZu • 

JiM\ Ji 

The total convective heat flux to the surface is given as: 

Qconv = ^ t ^ I't "*■ K v VT v + X h l pD^ 7 Xi • (5) 

Gas-Surface Interactions 

Two finite-rate, gas-surface interaction models are considered in the present 

study. 

Park Model 

Three kinds of gas-surface interactions for the carbonaceous ablating surface are 
considered. The gas-solid reactions include oxidation, Eqs.(6)-(7), and nitridation, 
Eq.(8). Following the assumption made in Park’s work, 7 ' 9 the surface catalysis reactions 
are assumed to be negligibly small for both oxygen atom and nitrogen atom 
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recombination at the ablating surface. The sublimation of carbon produces primarily C 3 
by Eq.(9). The species C 5 and C 7 are not included, because the entry conditions studied 
in this work have a stagnation point pressure of less than 1 atm. Thus, the carbon mass 
blowing rates due to gas-surface interaction are: 


m \ = pCo v oPo .f 

Mq 

(O + C( s ) —> CO) 

(6) 

M r 

m 2 =2pC 0 v 0 (3 0 c 

Mq 2 

(O 2 2C(_y) — ^ 2 CO) 

(7) 

m 3 = P c N v nPn , , 
M n 

(N + C (S) ^CN) 

(8) 

m 4 = p(C Cj E ~ C c 3 Wc 3 Pc 3 

QC(s) Q) 

(9) 

m c = mi + + m 4 - P( s )S 


(10) 


Here, v l is defined as^kT w !2mn l . The total carbon mass blowing rate, m c , and the 

carbon surface recession rate, S , can be computed from the total carbon mass loss. J3 is 
the efficiency of each surface reaction. 9 po is given as 0.63exp(-l 160/7V) . pn and P 02 are 
assumed to be 0.3 and 0.5, respectively, and p C3 is set to 1. The saturated carbon vapor 
pressure (Pg) for carbon-phenolic is expressed as 9 : 

( 90908 

p E =6.27xl0 15 e Tw 
and for carbon as: 

( 90845 ) 

p E =5.19xl0 15 e Tw . 


Species mass conservation at the surface is written as: 

- pD^Xi + Pv w C t = N t + m g C ig (1 1) 

The first term on the left-hand side is mass transfer through diffusion, and the second 
term is mass transfer due to convection. On the right-hand side are the source terms due 
to gas-surface interaction and pyrolysis gas injection. For graphite or pure carbon heat- 

shield materials, the pyrolysis gas injection rate, m g , is zero. Here the source term, N , , 

is defined as: 
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(12) 


N co ~ 


M 


CO 


Mr 


+ m 2 


M 


CO 


Mr 


> N CN ~ ™3 


M 


CN 


M c 


~ . M AT * . M f) * 

N N =~ m 3r—, N 0 =-m l ——, N 0 =-m 2 
M Q M Q 

and N, = 0 for i = C0 2 , N 2 , NO, C 2 , and C. 


, Nq = m 4 , 

Mo, 

2 M c ’ 


Based on global mass balance, the following equation for total mass blowing is 
expressed: 


pv w =m c +rh g (13) 

Zhluktov Model 

The following surface-kinetic model is suggested in the work of Zhluktov and Abe 10 : 

1. O + (C) <=> (C-O) 

2. 0 2 + 2(C) <=> 2(C-0) 

3. 0 2 + (C) <=> (C-O) + O 

4. C0 2 + (C) <=> (C-O) + CO 

5. (C-O) <==> CO + (C) 

6. O + (C-O) <==> C0 2 + (C) 

7. 2(C-0) <=> C0 2 + 2(C) 

8. (C) <==> C + (C) 

9. 2(C) <==> C 2 + 2(C) 

10. 3(C) <=> C 3 + 3(C) 

1 1 . N + (C) <=> (C-N) 

12. (C-N) + N <==> N 2 + (C) 

The corresponding reaction rates are as follows: 

r \ = k f\ (pZo®° ~®o /K 1 ) 
r 2 =k r 2 (K 2 pzo 2 (®°) 2 -(®o) 2 ) 
r 3= k f3 ( PXo 2 ®° ~ PXo ®0 1 K 3 ) 
r 4 = k f 4 ( PXC0 2 ®° ~ PXC0®0 1 K 4 ) 
r 5 =kf 5 (® 0 -PXC0®° /K 5) 
r 6 ~ kf 6 ( PX0®0 ~ PXC0 2 ©° 1 K 6 ) 

(14) 
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r i = k fiii®o ) 2 ~PZco 2 (©^) 2 / K'j) 

r 8 =^ r 8©°(-^8 ~ PXc ) 

r 9 =k r 9 (Q°) 2 (K 9 -pZc 2 ) 

TO = ^rl0(©°) 3 (-^10 ~PZC 3 ) 
r ll = kf n (pZ N ®° ~®N /^ll) 
r 12 = k rl2( K 12PXN®N -PXN 2 ® 0 ) 


It is assumed that two possibilities exist for K) and Kj j. The first is mobile adsorption: 


K 


L = B kT { 2mnikT lt 2 -T d . / T 

r n ' T ' 5 


atm 


and the second is immobile adsorption: 


J_ kT 2 rnnjkT 3 / 2 -T d . / T 

Kt V A 2 


atm 


(15) 


(16) 


where i = 1 or 1 1, P 0 = 1.01 325X1 0 5 Pa, and B = 3.5 X 10 19 m' 2 . 

The equilibrium constants Ki to are not independent. They are related via gas phase 
equilibrium constants of corresponding dissociation processes: 

K 2 = (KpKoj 
K 3 = K x Kq 2 
K 4 = K\K COl 
K 5 =K % I{K x K co ) 

K 6 =K z I{K x K co K CO i ) (17) 

K 7 = K $/(( K l ) 2 KcoKco 2 ) 

K 9 ={K Z ) 2 /Kc 2 

K 10 = (K s ) 2 /(Kc 2 Kc 3 ) 

K n =\/{K n K N2 ) 


The recommended expressions for necessary forward or backward reaction rate constants 
as follows: 


7 



k f\ = £ \ f o 


trl =e 2 B(^)e- T ‘i ,T 

h 

kf3=eiFo 2 e- T °>' T 

kf A = £ 4 F C0 2 

k f5 =s 5 B(~)e~ r ^ /T 

k f6 =e 6 F 0 e- T °6 /T 
k f7 =e 7 B(^)e- T ° 7 /r 


( 18 ) 


k rS ~ £ & F C 
k r 9 = £ 9 F C 2 

k rlO = £ \0 F C 3 
kf\\ = £ \\ F N 

k r 12 = £ n F N 2 e TalllT 
where 

Fi = Pq / ^ImrtikT 

£ \ =l,e 2 = 0.0008, = U 4 = 0.9, 

£■5 = 0.1, £g = 0.8, e 7 - l,£g = 0.24, 

£ 9 = 0.5,£ - io = 0.023, j = 1, fj2 =1, 

T dx =45,000, 

T d \ j = 36,600 

F a2 = 2T d] -T D q 2 
T al = F D0 2 ~ T d\ 

T a\2 =F DN 2 ~ F d\\ 

T a5 =40, 000, 34, 000, or 27, 000; 

T a6 = 2000 or 1000 

T al = 40, 000, 34,000, or\ 3, 700. 

The rates of species production on the surface are: 
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™0 =(- r l + '3 -reWo 
™C0 = ( r 4 + r 5 )M CO 
m C0 2 = (~ r 4 + r 6+ r 7 W CO-2 

™C =r $ M C 
m Cl =r 9 M Cl 

™C 3 = r lO m C 3 (19) 

m NO =0 
m C N =0 

™N =(-'11 ~ r \l) M N 
™0 2 = (~ r 2 ~rj,)Mo 2 

™N 2 = r \2 M N 2 

For a stationary regime we have: 

™(C-0) ~ ('ll + 2 '2 + r 3 + '4 ~ '5 _r 6 - 2^7 )M o -0 
m (C-N) ~ ( r l 1 ~ '12 Wn = 0 • 


Using these expressions, the total surface mass blowing rate is, 

m = Zfhi =M c (r 5 +r 6 +r 7 +r 8 +2 r 9 +3r 10 ). (21) 

The sum of surface coverage concentrations is equal to 1 ; that is, 

<=>0+©JV+©° =1 (22) 

Species conservation at the surface is written as: 

— pDfi ' Xi T P^w^i ~ m i + riigCj g (23) 


The total number of unknowns is 13, including: 

[zo,Z0 2 ,ZN,ZN 2 ,ZC,ZC 2 ,ZC l ,ZCN,ZC0,ZC0 2 ] 


The nonlinear equations can be solved through iterations by the following equation, 

I 

Fi 


AZj = -| 


( dF t ^ 1 


y dz jj 


(24) 
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The surface temperature (T w ), pyrolysis gas injection rate (m g ), and species 

concentrations of pyrolysis gas (C ig ) must be specified. The pyrolysis gas is assumed to 

be at chemical equilibrium before being injected into the main stream. The surface 
temperature and pyrolysis gas injection rate are obtained from the material response code 
through iterations. Species concentrations, C i w , and velocity, v w , at each surface grid 

point are updated in each sweep (both backward and forward) until a steady-state flow 
solution is reached. Thus, the carbon mass injection and surface recession rates can be 
obtained as part of the flow solution. 

Results and Discussion 

The general carbon-phenolic ablation model without pyrolysis gas injection can 
be used to simulate graphite ablation. In the first test case, the interaction between a 
graphite model and a typical arc-jet stream conducted in the Interactive Heating Facilities 
at NASA Ames Research Center is studied. The stream total enthalpy is estimated to be 
approximately 27 MJ/Kg, and the measured stagnation point pressure is 0.75 atm. The 
free stream velocity is 5354 m/s, the density is 0.003 kg/m 3 , and the temperature is 1428 
K. The free stream species mass concentrations are C 02 = 0.0, Cn 2 = 0.6169, C N o = 
0.0046, Cn = 0.1212, and Co = 0.2573. In this arc-jet stream, oxygen is fully dissociated 
and nitrogen is partially dissociated. The graphite model is a 10° half angle sphere-cone 
with nose radius (Rn) of 1 .905 cm. The geometiy and computational grid for this test 
model are shown in Fig. 1 . The surface temperature obtained from a solid-fluid coupled 
chemical equilibrium ablating surface computation is used as the temperature boundary 
condition in the finite-rate ablation surface calculation (Fig. 2). This is to study the 
surface reaction rates of Park and Zhluktov under the same given surface conditions 
without a fully coupled computation. Additionally, the predictions based on the steady- 
state ablation surface conditions will also be presented and compared in the final paper. 
The predicted carbon mass blowing rate distributions over the model surface are shown 
in Fig. 3. The comparison of predicted surface heat fluxes is presented in Fig. 4. 

The second test case simulates the flow over a Phenolic Impregnated Carbon 
Ablator (PICA) heatshield at estimated Stardust re-entry peak heating conditions. PICA 
has density of about 240 kg/m 3 , which is much lower than fully dense carbon-phenolic, 
which has a density of 1440 kg/m 3 . The aeroshell of the Stardust re-entry vehicle is a 60° 
half angle spherical cone with the nose radius (Rn ) of 0.2286 m, base radius (Rb ) of 
0.4064 m, and comer radius (Rc) of 0.02 m. The free stream velocity is 11137 m/s, the 
density is 0.000234 kg/m 3 , and the temperature is 238 K. The geometry and 
computational grid for the Stardust aeroshell are shown in Fig. 5. The surface 
temperature and pyrolysis gas injection rate distributions 11 used as part of the boundary 
conditions are shown in Fig. 6. The predicted total mass blowing rates, including char 
and pyrolysis gas, are presented in Fig. 7. Figure 8 shows the comparison of predicted 
convective heat fluxes. The stagnation streamline solution computed in Park’s study 8 is 
approximately 5 % higher than the current finite-rate prediction. 

The third test case predicts the flow over a fully dense carbon-phenolic heatshield 
at peak heating point of the proposed MSR EEV trajectory. The forebody heatshield is a 
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60° half angle spherical cone with a nose radius (Rn ) of 0.352 m, a base radius (Rb ) of 
0.45 m, and a comer radius (Re) of 0.02 m. The free stream velocity is 10100 m/s, the 
density is 0.0006634 kg/m 3 , and the temperature is 255 K. The peak heating stream 
velocity and the geometry of the MSR EEV aeroshell are similar to those of Stardust. 
MSR EEV has a larger Rb and Rn, and the stream density is almost 3 times as high as 
Stardust. The geometry and computational grid for this test case are shown in Fig. 9. The 
temperature and pyrolysis gas injection rate distributions used in this test case are shown 
in Fig. 10. The predicted total mass blowing rates (char and pyrolysis gas) are presented 
in Fig. 1 1 . In this test case, a large portion of surface mass blowing comes from the 
pyrolysis gas injection, which is approximately 0.09 kg/m 2 -s. The carbon char mass 
blowing rate is smaller than the pyrolysis gas injection rate over the entire heat-shield 
surface. Figure 12 shows the predicted surface heat fluxes. The calculated heating 
distribution using chemical equilibrium surface is quite different from using finite-rate 
surface. 

Further details of the predicted flow fields and the effect of the mass diffusion 
model on the ablation rate will be included in the final paper. 
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Figure 3 : Predicted carbon mass blowing rate distributions 
for case 1 
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Figure 4: Comparison of surface convective heat flux 
predictions for case 1 
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Figure 5: Geometry and computational grid for case 2 
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Figure 6: Surface temperature and pyrolysis gas 
blowing rate distributions used in case 2 simulation 
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Figure 7: Predicted total mass blowing rate distributions 
for case 2 



Figure 8: Comparison of convective heat flux predictions 
for case 2 
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Figure 1 1 : Predicted total mass blowing rate distributions 
for case 3 



Figure 12: Comparison of convective heat flux predictions 
for case 3. 
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